Spatial Prediction with INLA

Introduction

This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation.

GMRF Background

Blangiardo and Cameletti (2015) section 6.1.

• Observations aggregated to disjoint areal regions indexed by \(i\).
• Each region has unique parameter \(\theta_{i}\).
• \(\mathcal{N}(i)\) is the set of indices of neighbors of region \(i\) and \(\mathcal{N}_{i} = |\mathcal{N}(i)|\) is the number of neighbor of region \(i\).
• Local Markov property: given \(\boldsymbol{\theta}_{\mathcal{N}(i)}\), \(\theta_{i}\) is independent of all other \(\theta_{j}\).

• Then the precision matrix \(\mathbf{Q}\) of \(\boldsymbol{\theta}\) is sparse because only neighbors have nonzero coprecisions.

• Besag-York-Molliè model:
  • Exchangeable random effects \(u_{i}\) with intrinsic conditional autoregressive (iCAR) structure.
  • \(u_{i}|\mathbf{u}_{-i} \sim \mathrm{N}\left(\mu_{i} + \sum_{j} a_{ij} (u_{j} - \mu_{j}) / \mathcal{N}_{i}, \sigma_{u}^{2} / \mathcal{N}_{i}\right)\)
  • (What are the \(a_{ij}?\)).
  • iCAR is an improper prior because covariance matrix not positive definite but this is ok for random effects.

SPDE Background

Find Lindgren et. al. (2011) and fill in motivation

• Spatial process \(\xi(\mathbf{s})\).

• Solve \((\kappa^{2} - \Delta)^{\alpha / 2}(\tau \xi(\mathbf{s})) = \mathcal{W}(\mathbf{s})\).
  • \(\kappa\) is a range parameter.
  • \(\Delta\) is the Laplacian.
  • \(\alpha\) is a smoothness parameter.
  • \(\tau\) a precision parameter.
  • Exact and sationary solution: \(\xi(\mathbf{s})\) is a Gaussian field with Matèrn covariance function.

• Finite element approximation $$.

Geostatistics Example

Toy dataset from Blangiardo and Cameletti (2015).

# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))

# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))

# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
toy_fit <- inla(
  y ~ -1 + intercept + f(spatial_field, model = spde),
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
toy_fit$summary.fixed

toy_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

Bei Dataset and inlabru

Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).

# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

bei_corners <- vertices.owin(Window(bei))
bei_domain <- cbind(bei_corners$x, bei_corners$y)
bei_full_mesh <- inla.mesh.2d(cbind(bei$x, bei$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_full_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei, pch = '.', col = 'blue')

bei_full_spdf <- as.SpatialPoints.ppp(bei)
# CHECK PRIORS!
matern_full <- inla.spde2.pcmatern(bei_full_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = matern_full) + Intercept
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf)

lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
n_quads <- 10
botleft <- cbind(runif(n_quads, 0, 950), runif(n_quads, 0, 450))
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
    cbind(
      botleft[r, 1] + c(0, 0, 50, 50),
      botleft[r, 2] + c(0, 50, 50, 0)
    )
  )})
bei_win <- do.call(
  union.owin,
  apply(botleft, 1, function(x){return(
    owin(x[1] + c(0, 50), x[2] + c(0, 50))
  )})
)
bei_hole <- bei[complement.owin(bei_win)]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)

plot(bei_hole, main = 'Observed Subregion', pch = '.', cols = 'black')

bei_hole_mesh <- inla.mesh.2d(cbind(bei_hole$x, bei_hole$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_hole_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei_hole, pch = '.', col = 'blue')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
# CHECK PRIORS!
matern_hole <- inla.spde2.pcmatern(bei_hole_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = matern_hole) + Intercept
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf)

lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

bei_samp_mesh <- inla.mesh.2d(cbind(bei_samp$x, bei_samp$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_samp_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
plot(Window(bei_samp), border = 'blue', add = TRUE)
points(bei_samp, pch = '.', col = 'blue')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
# CHECK PRIORS!
matern_samp <- inla.spde2.pcmatern(bei_samp_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = matern_samp) + Intercept
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf)

lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Bei Dataset without inlabru

centers <- gridcenters(dilation(bei_window_full, 40), 40, 20)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
                     count = NA_integer_, area = NA_real_)

for(r in seq_len(nrow(bei_df))){
  bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
                         bei$x < bei_df$x[r] + dx &
                         bei$y >= bei_df$y[r] - dy &
                         bei$y < bei_df$y[r] + dy)
  bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(bei_df$count, nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')

# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(bei_full_mesh, as.matrix(bei_df[bei_df$area > 0, c('x', 'y')]))

# Initialize exponential covariance structure for SPDE.
full_spde <- inla.spde2.matern(mesh = bei_full_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = full_spde$n.spde)
stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = bei_full_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
bei_full_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = full_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
bei_full_fit$summary.fixed

bei_full_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_pred <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- bei_full_fit$summary.linear.predictor[index_pred, 'mean']
post_sd <- bei_full_fit$summary.linear.predictor[index_pred, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(t(matrix(post_mean, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior Mean of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

plot(im(t(matrix(post_sd, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior SD of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

beihole_df <- data.frame(x = centers$x, y = centers$y,
                         count = NA_integer_, area = NA_real_)

for(r in seq_len(nrow(beihole_df))){
  beihole_df$count[r] <- sum(bei_hole$x >= beihole_df$x[r] - dx &
                             bei_hole$x < beihole_df$x[r] + dx &
                             bei_hole$y >= beihole_df$y[r] - dy &
                             bei_hole$y < beihole_df$y[r] + dy)
  beihole_df$area[r] <- area(Window(bei_hole)[owin(c(beihole_df$x[r] - dx, beihole_df$x[r] + dx), c(beihole_df$y[r] - dy, beihole_df$y[r] + dy))])
}

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beihole_df$count, nrow = length(unique(beihole_df$x)))), unique(beihole_df$x), unique(beihole_df$y), unitname = 'meters'), ncolcours = range(beihole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei, pch = '.', col = '#00000040')

References

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.

Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.

Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.